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We compare simulations of the Lyman-a forest performed with two different hydro- 
dynamical codes, Gadget-2 and Enzo. A comparison of the dark matter power 
spectrum for simulations run with identical initial conditions show differences of 1-3% 
, at the scales relevant for quantitative studies of the Lyman-a forest. This allows a 

^ ' meaningful comparison of the effect of the different implementations of the hydro- 

OO I dynamic part of the two codes. Using the same cooling and heating algorithm in 

CO i both codes the differences in the temperature and the density probability distribution 

' function are of the order of 10 %. The differences are comparable to the effects of 

boxsize and resolution on these statistics. When self-converged results for each code 
I are taken into account the differences in the flux power spectrum ~ the statistics most 

widely used for estimating the matter power spectrum and cosmological parameters 
from Lyman-a forest data - are about 5%. This is again comparable to the effects of 

• boxsize and resolution. Numerical uncertainties due to a particular implementation of 
Oh| solving the hydrodynamic or gravitational equations appear therefore to contribute 

only moderately to the error budget in estimates of the flux power spectrum from nu- 

• merical simulations. We further find that the differences in the flux power spectrum for 
^ I Enzo simulations run with and without adaptive mesh refinement are also of order 5% 

. or smaller. The latter require 10 times less CPU time making the CPU time require- 

ment similar to that of a version of Gadget-2 that is optimised for Lyman-a forest 
simulations. 
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1 INTRODUCTION 

There is now a well established paradigm for the origin of 
the Lyman-Q forest, the ubiquitous absorption lines due to 
neutral hydrogen in the spectra of high-redshift quasars. 
The absorption blue-wards of 1216 A is predominantly due 
to density fluctuations in the intervening warm (~ 10* K) 
photoionized inter-galactic medium (I GM) on scales larger 
than the Jeans length of the gas (see iRauchI il998t) for a 
review). Numerical simulations were inst rumental in estab - 
lishing the new pa r adigm in the 1990s dCen et alJ il993). 
IZhang et all Jl995l). iHernQuist etU] (Il99dl . iTheuns et alJ 
Jl998f) . IZhanget all Jl997t) '). The Lvman-g forest and the 
associated metal absorption probe the thermal and ioniza- 
tion history of the IGM as well as the interplay of galaxies 
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and the IGM from which they are formed. More recently 
the Lyman-a forest has also been established as a means of 
quantitative measurement of the underlying matter distri- 
bution and t hus a variety of cosmolog ic al paramet e rs (e.g . 
Croft et al.l lll998ll ICroft et alj J2002I) IViel et al.l J2004ii 
Sebak et alJ ^2005^ . IViel fc HaehneITll200eD . ISeliak et alJ 
J200d) : Vie". Haehnelt & Lewis (2006)). Numerical simula- 
tions thereby play a crucial role in inferring the linear mat- 
ter power spectrum and other derived quantities from the 
Lyman-Q forest data. With increasing sample sizes statisti- 
cal errors of measurements of the flux distributions have 
reached the percent level and the error budget is domi- 
nated b y systematic uncertainti es (Viel, Haehnelt & Springel 
(2004). McDonald et all (1200511 ). Uncertainties due to nu- 
merical simulations contribute signiflcantly to the error bud- 
get and the accuracy with which the flux distribution for 
given input physics and cosmological parameters can be 
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simulated has become important. Most studies so far have 
used convergence tests to assess uncertainties due to the nu- 
merical simulations and direct comparisons of cosmological 
hydrodynamical simulation performed with different codes 
have been rare. The differences between hydro-dynamical 
simulations of galaxy clusters with a wide range of different 
codes/methods have been studied in th e Santa Barbara clus - 
ter project (Frenk et al. 1999). Recentlv lO'Shea et alJ feOOST) 
performed a comparison between the grid based adaptive 
mesh refinement (AMR) code Enzo^ and the smoothed par- 
ticle hydrodynamics (SPH) code Gadget-2^. However, lit- 
tle has been done in this respect for h ydrodynamica l simu - 
lations of the Lyman-a forest data (see lTheuns et al] (Il998f) 
for a notable exception of a comparison between two SPH 
codes) Some comparisons of hydrodynamical simulations 
with approximate simulations of the Lyman-a forest data 
have been performed by McDonald et al. (20051 IZhan et alj 
^2005^ and Viel, Haehneh & Springel (2006). We present 
here a comparison of hydrodynamical simulations of the 
Lyman-a forest with Enzo and Gadget-2 which concen- 
trates on the statistical properties of the flux distribution. 
We are therefore mostly interested in properties of the mod- 
erate to low over-density gas which is responsible for Lyman- 
Q forest absorption. Of particular interest is the probability 
distribution of the gas density, the temperature, the result- 
ing flux distribution and the flux power spectrum. A major 
difference between grid-based and SPH codes is their treat- 
ment of shocks and their effects on the temperature distri- 
bution. We will also examine these differences. 

The plan of the paper is as follows. In describe 
the Enzo code and the Gadget-2 code and the different 
ways in which the codes solve the gravitational and hydro- 
dynamics equations. In H2.4I we describe the simulation set 
used in the comparisons. In i|3|we investigate how physical 
properties of both codes compare, in particular the gas dis- 
tribution and the 1-D flux power spectrum. Finally in tl3.3l 
we will look at the performance of each code in terms of 
CPU time consumption. 



2 HYDRODYNAMICAL SIMULATIONS OF 
THE Lyman-Q FOREST 

2.1 Grid-based simulations vs. SPH simulations 

The physical state of the gas responsible for the Lyman- 
a forest is largely governed by the competing processes of 
photoionisation and adiabatic cooling due to expansion. The 
low over-de nsity gas obeys a t ight temperature - de nsity re- 
lation, Ce.g. lKatz et al.llll996l) . lHui fc GnedinI lll997ll ) which 
can be approximated by 

T.r.(i)- a, 

where pb and P(, are the baryon density and mean baryon 
density respectively, T is the temperature. To and 1 < 7 < 
1.6 are parameters which depend on the reionisation his- 
tory model and the spectral shape of the UV background. 
Typical temperatures of the photoionized IGM are in the 
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^ http:/ /www. mpa-garching.mpg.de/gadget/ 



range 10000-20000 K. The optical depth for Lyman-a ab- 
sorption is proportional to the neutral hydrogen density, 
jGunn fc PetersonllTge^) . which, since the gas is in pho- 
toionisation equilibrium, is proportional to the square of the 
density times the recombination rate, 

-if:) 

where p — 2.7 — O.77. The factor A depends on the redshift, 
baryon density, temperature at the mean density, Hubble 
constant and the photoionisation rate. The optical depth is 
a faithful tracer of the matter distribution on scales larger 
than the Jeans length of the photoionized IGM. Even though 
the density fleld is only mildly non-linear on the relevant 
scales the thermal effects, the peculiar velocities and the 
non-linear relation between flux and optical depth make hy- 
drodynamical simulations mandatory for accurate predic- 
tions of the statistical properties of the flux distribution in 
the Lyman-a forest. Cosmological hydrodynamical simula- 
tions come in two basic flavours, SPH (Lucy 1977, Gingold 
& Monaghan 1977) and grid based codes. SPH codes like 
Gadget-2 use particles to represent the baryonic fluid. SPH 
is a Lagrangian method and hence the resolution is concen- 
trated in regions of high density. Grid based codes like Enzo, 
use a grid of cells to represent the gas properties. One of the 
options in Enzo is an adaptive reflnement of the grid where 
the grid resolution is increased in regions of high density re- 
sulting in a large dynamic range. Both codes have been used 
to stud y the Lyman - a fore s t (e.g Viel, Haehnel t fc Springel 
(2004) , llTvtler et alj J2004l) . iBolton et alj J2005l) . I Jena et all 
boOSi) . iMcDonald et all tOO^ ). SPH simulations of the 
Lyman-a forest in particular have been very successful but 
in principle one may think that a grid based code could 
offer better resolution of the low density IGM than SPH 
codes because in SPH simulations the resolution in low den- 
sity regions decreases when the gravitational clustering of 
the matter distribution becomes non-linear and high den- 
sity regions start to form. Grid-based codes could also offer 
a more accurate treatment of shocks which are relevant for 
the thermal state of the IGM. We will now briefly describe 
the methods implemented in Enzo and Gadget-2 to solve 
the gravitational and hydrodynamical equations. 

2.2 Enzo 

Enzo is a Eulerian adaptive mesh refinement code origi- 
nally developed by Greg Bryan and Mike Norman at the 
University of Illinois (Bryan fc Norman 1995b, Bryan fc 
Norman 1997, Norman fc Bryan 1999, O'Shea et al. 2004). 
The hydrodynamics solver employs the Piecewise Parabolic 
method combined with a non-linear Riemann solver for 
shock capturing. The Euleria n AM R scheme was first de- 
veloped by Bcrg cr fc Qligciji 1^8^ and later refined by 
iBereer fc ColeUal lil989l) to solve the hydrodynamical equa- 
tions for an ideal gas. Bryan fc Norman 1997 adopted such 
a scheme for cosmological simulations. The gravity solver in 
Enzo uses a N-Body particle mesh technique (Efstathiou et 
al. 1985, Hockney fc Eastwood 1988). 

We have used the publicly available version of 
Enzo (Enzo -1.0.1) which we have modified to use the 
Gadget-2 equilibrium chemistry solver as discussed in 

mm 
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Table 1. Spatial parameters and spatial dynamic range (SDR) for 
the Enzo AMR simulations, the Enzo static grid simulations and 
the Gadget-2 simulations with and without star formation. The 
particle numbers for Gadget-2 and the rootgridsizc for Enzo is 
shown in the loft hand column. In the Gadget-2 simulations the 
number of DM particles and SPH particles is equal and for the 
Enzo simulations the number of DM particles is equal to the 
rootgridsize. 



2.2.1 Enzo gravity solver 

The gravity solver in Enzo employs an adaptive particle- 
mesh algorithm. The potential is solved on a periodic root 
grid using Fast Fourier Transforms. In order to accurately 
account for the sub grids a multi-grid relaxation technique is 
used (see e.g. Norman & Bryan 1999). The force resolution 
is typically twice as coarse as the grid resolution. 



2.2.2 Enzo hydrodynamics solver 

Enzo uses the Piecewis e parabolic method (PPM), 
JWoodward fc Colellal Il983) for solving the hydrodynamic 
equations. A complete description of this method is not 
possible her e and we will only give short description (see 
iBrvan et alJ (Il995ll for more details) . PPM is a higher order 
accurate version of Godunov's method with a third order 
accurate piecewise parabolic monotonic interpolation and a 
non-linear Riemann solver for shock capturing. The method 
is second order accurate in space and time and explicitly 
conserves energy mass flux and momentum. It uses a dual 
energy formalism which allows the calculation of both the 
thermal energy and the total energy of the gas at each time 
step. This ensures a correct internal energy and the correct 
entropy jump at shock fronts and the correct temperature 
and pressure in hypersonic flows. This represents a major 
difference with respect to SPH codes which employ an arti- 
ficial v iscosity to capture shocks (see ISpringel fc Hernauisli 
J2002D '). 

In addition to solving the ideal gas dynamics Enzo also 
has several cooling and heating routines. For the cooling 
both non-equilibrium and equilibrium cooling functions are 
available. 

We will use here a modified version of Enzo, which. 



uses the Gadget-2 equilibrium chemistry solver for hydro- 
gen and helium with a uniform ultra-violet (UV) background 
based on the models of lHaardt fc Madaul ilQQST) . 



2.2.3 The adaptive mesh refinement 

The AMR ability of Enzo introduces finer and finer grids 
into areas of high density thus allowing maximum resolution 
where it is actually needed at a minimum computational 
expense. This ability to dynamically refine the resolution 
is essential for accurately tracking the non-linear collapse 
of rapidly evolving density fields. For Lyman-a forest sim- 
ulations the resolution in high density regions is less im- 
portant and we investigate here simulations with and with- 
out AMR. With AMR the hydrodynamical equations are 
initially solved on a uniform grid, the solutions are moni- 
tored and the patches of the initially uniform grid are re- 
fined if certain refinement criteria are met. Parent grids 
then produce child grids. For cosmological studies a refine- 
ment factor of 2 is normally used. This means that a child 
grid will have cells which have twice the spatial resolution 
of the parent grid. It is also worth noting that grids at 
the same level will have the same timestep but that this 
timestep may be different for grids at a different level of re- 
finement. For Lyman-a forest simulations it is not obvious 
how important the use of the AMR option is. Most of the 
absorption especially at high redshift is by gas of moder- 
ate (over) density. Note that previous studies of the Lyman- 
a forest with grid-based codes have generally not used AMR 
methods (e.g. lJena et aj2qq5l.rMcDonald et aj2005f) . How- 
ever as shown ]bv l^e^^^T(l2004l) the few strong absorption 
systems caused by dense regions contribute significantly to 
the flux power spectrum at all scales. We will investigate 
this further in >I3.2.4I 



2.3 Gadget-2 

Gadget-2 (Springel 2005), the updated version of Gadget- 

1 (Springel 2001) is a parallel TreePM-SPH code. On the 
scales relevant for the Lyman-a forest Gadget-2 in its 
TREEPM mode is similar to a PM code with some extra 
resolution due to the Tree part of the algorithm on small 
and intermediate scales. The gravitational components of 
Gadget-2 and Enzo are therefore somewhat similar on 
large scales most relevant for the Lyman-a forest. The hy- 
drodynamical components are, however, very different. 

We have used a version of Gadget-2 which is similar 
to the publicly available version as of August 2006. The only 
exceptions are that we have used a Gadget-2 equilibrium 
cooling algorithm supplied to us by the authors of Gadget- 

2 , and that some of the simulations discussed in sections 
t|2.4l and H3.2.4I were run with a version of Gadget-2 op- 
timised for speed for Lyman-a forest simulations where gas 
with an overdensity > 1000 and a tempera ture < 10^ K is 
turned into coUisionless star particles (see IViel et alj|2004l 
for more details). 



2.3.1 The Gadget-2 gravity solver 

We have used a version of Gadget-2 which employs a 
TreePM algorithm to solve the gravitational equations 
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Figure 1. Left Panel: Top: The DM power spectrum of the Enzo AMR simulations with a rootgridsizo of 200^. The arrows indicate 
the Nyquist frequency for the rootgrid for the 60 Mpc , 30 Mpc and 15 Mpc boxes (comoving) from left to right. Bottom: 
The fractional difference between the Enzo simulations with and without AMR [(Enzo(AMR) - Enzo(static) ) / Enzo(AMR) ]. 
Right Panel: The fractional difference between the DM power spectrum of simulations with Gadget-2 and Enzo(AMR) [(Gadget-2 - 
Enzo{AMR) ) / Gadget-2 ]. All results are for z = 3.0. 



(Xu 1995, Bode et al. 2000, Bagla & R ay 2003). The 
Tree PM algorithm is a hybrid of the tree jBames fcJIud 
[l98d) and particle mesh m ethods iEfstathio^^^Lr^SSl 
iHockney fc EastwoodllT988r) . It utilises the best elements of 
both making the gravitational force determination more ac- 
curate and efficient. The potential is split into two compo- 
nents <j)k = (/>)°°'^ + <?if where 

"^1°"*^ = (t>kexp{-k'^r1) (3) 

and Ta is the spatial scale of the force spHt which is 
usually set to a little larger than the mesh spacing. For the 
simulations performed here we used the Gadget-2 default 
value of rs equal to 1.25 times the mesh spacing and a mesh 
spacing equal to the cube root of the total number of dark 
matter particles. The long range force is then computed us- 
ing mesh methods making the long range force almost exact. 
The short range force is computed using a tree algorithm 
which calculates the gravitational particle-particle forces on 
small and intermediate scales in an efficient manner. 

In order to conserve the symplectic nature of the leap- 
frog time integration for the case of individual timesteps 
the Hamiltonian is separated into a kinetic part and a long 
range and short range potential. Gadget-2 then evolves all 



particles using individual timesteps hence reducing the com- 
putational overhead that would be associated with evolving 
all the particles using the minimum allowed timestep. The 
splitting of the time integration is similar to what is done in 
the TreePM algorithm (see^pringcl 2005 for more details). 
TreePM codes offer an excellent compromise between speed 
and accuracy. 



2.3.2 The Gadget-2 Hydrodynamics solver 

The Gadget-2 hydrodynamics solver uses the SPH formal- 
ism. SPH can be thought of as a discretisation of a fluid 
which is then represented by particles. Continuous fluid 
properties are then defined using kernel interpolation. The 
particles sample the gas in a Lagrangian sense thus making 
SPH methods very powerful for following structure forma- 
tion in cosmological simulations. 

The thermodynamic state of each fluid element can ei- 
ther be defined in terms of its thermal energy per unit mass, 
Ui, or its entropy per unit mass, Si. In Gadget-2 the entropy 
per unit mass, Si, is used (see Springel & Hernquist(2002)). 
The code conserves both energy and entropy even when fully 
adaptive smoothing lengths are used. This represents a ma- 
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Figure 2. Left Panel: A comparison of the Enzo(AMR only) and Gadget-2 gas temperature over-density relation at z = 3 for the 
(15, 200) simulations. The diamonds represent the Enzo simulation while pluses represent the Gadget-2 simulation. 10000 values from 
each simulation are plotted. Right Panel: PDF of the volume-weighted gas temperature distribution for the (15, 200) simulations with 
Enzo with and without AMR and Gadget-2 with and without simplified star formation. All the results are a,t z = 3.0. 



jor change in methods b etween Gadget- 1 a nd Gadget- 
2 which is investigated in lO'Shea et all ^0^. 

A potential drawback of SPH is the approximate way 
in which it c a pture s shocks by use of an artificial viscosity 
fsee lSpringell (|200^ for more details). 

The large differences in the methodology of the hydro- 
solvers make a comparison of Lyman-a forest simulation 
with both codes very interesting for a test of the sensitivity 
of these simulations to a particular numerical method. 



2.4 Simulation parameters 

We have performed simulations with parameters of the 
concordance cosmological model with = 0.74, fim = 
0.26, fib = 0.04 63,(78 = 0.85,71 = 0.95 and h = 0.72. 
T hese simula t ion pa rameters correspond to the 'B2' model 
of lViel et all 112004). Identical initial conditions were used 
for the simulations with both codes. Enzo was run with the 
implementation of the Gadget-2 equilibrium chemistry so 
that cooling a n d hea ti ng were also treate d in the same way. 
iTheuns et alJ (Il998l). IViel et alJ i2004l). iMcDonald et'al] 
j2005h . lJena et al.l (|2005^ and lBolton et all ^2005^ have Tii^ 
vestigated the effect of box size and resolution of hydrody- 
namical simulations on a variety of statistics of the Lyman- 
Q forest flux distribution most importantly the flux proba- 
bility distribution, the effective optical depth and the flux 
power spectrum. Unfortunately it is currently not possible 
to run hydrodynamical simulations for which these statis- 
tics are fully converged especially not for a large parameter 
space. Convergence and box size studies are therefore essen- 
tial for quantitative studies of the Lyman-a forest. Gener- 
ally compromises have to be made and the application of 
corrections and an analysis of the corresponding errors are 



necessary. Our aim is here to investigate how the uncer- 
tainties between codes employing different hydrodynamical 
methods compare to other errors in a quantitative analysis. 
We have thus ran simulations with up to four different num- 
bers of basic resolution elements (number of particles and 
grid cells, respectively) and for three different box sizes. The 
Enzo simulations were run without AMR and with up to 4 
levels of grid refinement. We have chosen a refinement level 
of 4 so as to match Spatial Dynamic Range (SDR) of SPH 
calculations currently used in calculating the fiux statistics 
of the Lyman-a forest. The simulation parameters are sum- 
marised in Table Q Note that only the higher resolution 
simulations resolve the Jeans mass well. We will come back 
to this later. 

The SPH nature of the Gadget-2 simulations leads to 
a varying resolution similar to that of an AMR code. In 
order to get a feel for how the resolution of the different 
simulations compare, the last column of Table gives the 
SDR. For the SPH simulation the SDR is calculated as SDR 
= L/e where e is the gravitational softening length and is 
calculated before the simulation begins by dividing the box- 
size L by the number of particles along one axis times some 
constant factor. 

The Enzo simulations were run with a static grid and 
with AMR. For Enzo the SDR is calculated as SDR = 
A^root X 2', where A^root, is the size of the root grid in ID 
and I is the refinement factor. In the static grid simulations 
the grids are fixed throughout the simulation without any 
refinement. As a result the spatial resolution in high den- 
sity regions will be comparatively poor in these simulations. 
Most of the absorption in the Lyman-a forest is, however, 
produced by regions of low or moderate density. As we will 
see later the differences in the statistics of the flux distribu- 
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Figure 3. PDF of the volume- weighted gas density distribution for simulations with Enzo with and without AMR and Gadget-2 with 
and without simplified star formation. The panels are for simulations with a boxsizc of 15/i^^,30/i^^ and QOh^^ Mpc (comoving) from 
left to right. All the results are a.t z = 3.0. 



tion between simulations with and without AMR are there- 
fore actually moderate. The static grid simulations would be 
comparable in resolution to a Gadget-2 simulation with a 
softening length equal to the mean inter-particle spacing. 

Unfortunately the improved resolution of the AMR sim- 
ulations comes at the expense of a significant increase in 
computational time. For the AMR simulations we set the 
maximum refinement level to 4 beyond which an artificial 
pressure support is introduced to prevent further collapse. 
We thereby experimented with the values of the parame- 
ters for the minimum pressure support and checked that 
the thermodyamic properties of the cells in question were 
not affected. The mesh-refinement criterion was set to the 
standard values of 4 for the dark m atter (DM) and 8 for 
the baryons (see lO' Shea et al"! l|200^ for more details). This 
means that a grid will refine when its DM density reaches 
a factor 4 greater than the root DM density or when its 
baryonic density reaches a factor of 8 greater than the root 
baryonic density. To make contact with the simulations used 
in actual measurements of the matter power spectrum from 
Lyman-a forest data we also inves tigated some of the sim- 
ulations used in IViel et alJ l|200^. These simulations em- 
ploy a simplified star formation criterion that turns all gas 
with an over-density, with respect to the mean baryonic den- 
sity, > 1000 and temperature < 10^ K into colhsionless star 
particles. This substantially reduces the required computa- 
tional time by eliminating the short dynamic time scales 
associated with high density gas and significantly speeds up 



Lyman-a forest simulations with Gadget-2. The effect on 
the statis tics of the fiux di stribution has been shown to be 
small (see lViel et al] ^M^j) and we have labelled these sim- 
ulations as Gadget-2 (stars). 

Most of these simulations have very large particle numbers 
(2 X 400"^). Unfortunately for these simulations only a com- 
parison with the Enzo static grid simulations is feasible with 
our limited computational resources. We will concentrate 
our comparison on simulation outputs at z = 3 the centre 
of the redshift range 2 < z < 4 relevant for quantitative 
measurements of the matter power spectrum studies from 
Lyman-a forest data, but will briefly discuss simulations of 
the Lyman-a flux power spectrum at z = 2 and z = 4 in 

mm 

3 CODE COMPARISONS 

3.1 The dark matter distribution 

We start with a comparison of the DM distribution. O'Shea 
et al. (2005) have recently performed such a comparison 
and found moderate differences but note that we are in- 
terested in a different application of the code than investi- 
gated in O'Shea et al. (2005). The most relevant statistical 
property of the matter distribution is the power spectrum 
P(k) =< >, where is the Fourier transform of the 

density field. O'Shea found very good agreement at large 
scales with deviations at small scales. 
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Figure 4. PDF of the flux distribution for simulations with Enzo with and without AMR and Gadget-2 with and without simplified 
star formation. The panels are for simulations with a boxsizc of 15h~^ ,30h~^ and 60h~^ Mpc (comoving) from left to right with 
200^ particles/rootgridsize. The bottom panel shows the fractional differences between the Enzo(AMR) simulation and the Gadget- 
2 simulation (solid curve), the Enzo(AMR) and Enzo(static) simulations(dotted curve) and between the Gadget-2 simulations with 
and without star formation (dashed curve in the left panel) All the results are at 2 = 3.0. 



The left panel of Figure shows the DM power spec- 
trum of the Enzo(AMR) simulation in the form A^(fc) = 
P{k)k'^ for three simulations with a 200"^ root grid but dif- 
ferent box sizes at z = 3. The bottom panel in the left panel 
of Figure0shows the percentage difference between simula- 
tions with and without AMR in the form [(Enzo(AMR) - 
Enzo(static) ) / Enzo(AMR) ]. As expected the AMR 
simulations show more small scale power as the particle 
mesh algorithm becomes more accurate at small scales due 
to mesh refinement. The differences at large scales are very 
small, of the order of 1%. In the right panel of Figure Q we 
have plotted the fractional difference between the DM power 
spectrum of the Enzo(AMR) and the Gadget-2 simula- 
tions in the form [(Gadget-2 - Enzo(AMR) )/Gadget- 
2 ]. Similar to O'Shea we find differences of ~ 1% at large 
scales. 

The strong increase at small scales is due to the some- 
what different reso l ution limits of the simulations com- 
pared. |o]She^^^S] ll2005l) came to a similar conclusion and 
showed that Enzo and Gadget-2 produce very similar re- 
sults, even at small scales, when a low over-density thresh- 
ol d is used for the mes h- refinement criteria. As pointed out 
bv lO'Shea et al.l (|200^, and verified by us, Gadget-2 has a 
higher force resolution and hence more power at small scales 
due to the better force resolution of the Tree algorithm when 
simulations of similar SDR are compared and standard pa- 



rameters are used in both simulations. Note, however, that 
for simulations of the Lyman-a forest the differences at the 
relevant scales are very small. This makes a meaningful com- 
parison of the effect of the different hydrodynamics solvers 
on the statistics of the flux distribution - the main aim of 
this paper - possible. 

3.2 Properties of the gas distribution 

3.2.1 Shock heating and the thermal state of the gas 

As we described in at moderate to low over densities 
^ 5) the temperature - density relation is well approx- 
imated by a power law. In the left panel of Figure |5| we 
have plotted the temperature - density relation for a simu- 
lation with 2 x200^ DM and gas particles/rootgridsize in a 
15/i~^ Mpc box, hereafter labelled as a (15, 200) simulation, 
for Enzo(AMR) and Gadget-2. 10000 points are plotted 
for each simulation. This gives a feel for the level of shock 
heating produced by each code and also emphasises that the 
vast majority of the gas lies very close to a line representing 
the power law approximation. The temperature is volume- 
weighted in both cases. We have chosen volume-weighted 
temperatures as the flux statistics of the Lyman-a forest 
are volume-weighted. Temperatures and densities were cal- 
culated in the same way as for the mock absorption spectra. 
Note that the differences for the mass-weighted tempera- 
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tures were somewhat smaller. 

Overall the agreement between the two codes is remark- 
able given the very different ways in which both codes treat 
shocks. As mentioned previously Gadget-2 uses a conser- 
vative entropy formalism to treat shocks while Enzo uses a 
non-linear Riemann solver. In principle since Gadget-2 em- 
ploys an artificial viscosity to capture shocks Enzo should 
resolve shocks more accurately and one may have expected 
that weak shocks may occur in low density regions which 
Gadget-2 does not capture properly. This appears not to 
be the case. The amount of shock heated gas and its tem- 
perature distribution is very similar. This can be seen more 
clearly in the right panel of Figure |5| where we have plot- 
ted the volume weighted probability distribution function 
(PDF) of the temperature for simulations with Enzo with 
and without AMR and for Gadget-2 with and without the 
simplified star formation criterion. The differences are of 
order 10% and can be at least partially attributed to differ- 
ences in the PDF of the density which we will discuss in the 
next section. 



3.2.2 The probability distribution function of the gas 
density 

Figure 13 shows the volume- weighted PDF of the gas distri- 
bution for simulations with a range of box sizes. The agree- 
ment between the Enzo simulations with and without AMR 
in this linear volume-weighted plot which emphasises gas 
around the mean density is very good. The differences be- 
tween the simulations with Gadget-2 and Enzo are some- 
what larger, of order 10%. Note, however, that this is smaller 
than the differences due to changes in box size and resolu- 
tion. Overall the agreement is again very good. 

As demonstrated in the previous sections, a SPH code 
and a grid-based code differ in their resolution properties 
and it is not trivial to run simulations with the "same res- 
olution" because of differences in force resolution and the 
way the resolution is distributed between regions of differ- 
ent densities. The (small) differences between simulations 
with Enzo and Gadget-2 are thus not surprising. 



3.2.3 The probability distribution of the flux 

We have computed the flux distribution for 1000 random 
lines of sight through the simulation box. The optical depth 
has been rescaled in the standard way to match the observed 
effecti ve optical depth at z = 3 as given by ISchave et alJ 
i2003h . Tefi — 0.363. Figure |1] shows the corresponding 
probability distribution of the flux for simulations with 
Enzo with and without AMR and Gadget-2 with and 
without simplified star formation. Overall the flux distri- 
butions are very similar. Typical differences between the 
simulations with Enzo and Gadget-2 are 5-10%. The dif- 
ferences are again smaller than those due to changes of box 
size and resolution. The differences between the Gadget- 
2 simulations with and without star formation are less than 
3%. 
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Figure 5. The top panel shows the flux power spectrum for sim- 
ulations with Enzo(AMR) with a boxsize of 15ft~l,30/i-i,60/i~^ 
Mpc comoving. The middle panel shows fractional difEer- 
ences between the Enzo simulation with and without AMR 
[Enzo(AMR) - Enzo(static) )/Enzo(AMR) ]. The top and mid- 
dle panel are for 2 = 3 and the line styles are identical. The 
bottom panel shows the effect of resolution on the flux power 
spectrum at 2: = 2. The diff'erences are relative to the highest 
resolution run e.g. ((15, 400) - (15,200)) / (15,400). 



3.2.4 The flux power spectrum 

In Figure |S] we show the flux power spectrum for the 
Enzo AMR simulations for different boxsizes at z = 3. 
The middle panel shows the fractional difference between 
Enzo simulations with and without AMR. The solid curve is 
for simulations with a boxsize 15 h~^ Mpc, the dashed curve 
is for a boxsize of 30 h^^ Mpc box and the dotted curve is for 
a box size of 60 h'^^ Mpc. At large scales the differences are 
less than 4%. At the resolution limit the differences increase 
as expected. The effect of the AMR is most signiflcant at 
small scales. Note that the force resolution limit (twice the 
cell length) is about an order of magnitude off the graph on 
the right hand side. The bottom panel of Figure |K| demon- 
strates the convergence of the flux power spectrum by com- 
paring the flux power spectrum for the static grid simula- 
tions with a box size of 15 h~^ Mpc at z = 2. The differences 
between the (15,200) and (15,400) simulations are less than 
2 percent on the relevant scales suggesting that a resolution 
of 150/i~^ kpc (comoving) is required to reach co nvergence. 
This is in good agreement with the results by IViel et all 
l|20^ for the Gadget-2 simulations also used here. Note, 
however, that a similar comparison by Jena et al. (2005) 
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Figure 6. The fractional difference of the flux power spectrum 
of simulations with Enzo(AMR) and Gadget-2 for different box 
sizes [(Gadget-2 - Enzo(AMR))/Gadget-2 ]. Also shown is the 
difference for Gadget-2 simulations with and without star for- 
mation. The shaded region indicates the range of wave numbers 
used by Viel et al. (2004) to infer the linear dark matter power 
spectrum. 
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Figure 7. The fractional difference of the flux power spectrum 
of simulations with Enzo (static) and Gadget-2 for 400'' gas 
particle/rootgridsize simulations with different box sizes. The dif- 
ference is of the form [(Gadget-2 - Enzo(static))/Gadget-2 ]. 
All Gadget simulations are in this case run with the simplified 
star formation algorithm. The shaded region indicates the range 
of wave numbers used by Viel et al. (2004) to infer the linear dark 
matter power spectrum. 



for static grid Enzo simulations with the same resolution as 
shown in the bottom panel (their Figure 7) showed signif- 
icantly larger differences. A discrepant result for which we 
do not have an explanation. 

We have also investigated how the level of AMR refine- 
ment effects the Lyman-a flux power spectrum. The results 
for simulations with refinement level 2 lie in between those 
with refinement level 4 shown in Figure |3 and the static 
grid simulations. This suggests that a relatively high refine- 
ment level may be needed to correctly account for the effect 
of strong absorption systems on t he flux power spe ctrum 
caused by high density gas which IViel et al.l ^QM) have 
shown to extend to large scales. 

fn Figure|S|we show the fractional difference of the flux 
power spectrum of the Gadget-2 and Enzo simulations 
for the 200^ simulations. Apart from the smallest scales in 
the lowest resolution simulation the differences are about 
5%. The differences are scale dependent and appear to 
decrease with increased resolution. Unfortunately we did 
not have the computational resources available to run 
Enzo(AMR) simulation with a root-grid size of 400^. In 
Figure [7| we therefore show the differences between 400^ 
Gadget-2 and Enzo(static) simulations. The differences 
are again about 5%. The reader should thereby keep in mind 
the differences of up to 4% between the Enzo(AMR) and 
Enzo(static) simulations. Note that Viel et al. (2004) 
found the (60,400) simulations to be the best compromise 
between box size and resolution in their measurement of 
the matter power spectrum from Lyman-a forest data. 
We have also looked into the differences at z = 2 and 
z — 4 and found the differences to be redshift dependent. 
At z = 4 the differences are similar than those at z = 3 



while at 2 = 2 they are somewhat larger. By investigating 
the z — 3 simulations with different effective optical 
depth we verifled that the change of the differences of the 
simulated flux power spectra with redshift is mainly but 
not only due to the strong evolution of the effective optical 
depth. Also worth noting is the excellent agreement be- 
tween the Gadget-2 simulation with (thin solid curve) and 
without star formation (thick solid curve) shown in Figure|n] 



3.3 CPU time requirements 

We have performed a series of timing tests for both 
Enzo and Gadget-2 for a selection of the simulations used 
in this study. For the timings, simulations were carried out 
on the distributed memory Sun Cluster at the Institute of 
Astronomy in Cambridge and COSMOS a shared memory 
machine located at the Department of Applied Mathemat- 
ics and Theoretical Physics (DAMTP) in Cambridge. On 
the distributed memory machine the timing tests were per- 
formed on 32 single node 500 Mhz processors with 2 Gbytes 
of memory and a 100 Mbit ethernet connect. The latency 
of the connection is approximately 300 microseconds. Al- 
though these processors are relatively slow by today's stan- 
dards they should still provide an illustration of the relative 
performance of the codes in different configurations. The 
tests were performed for simulations with a boxsize of 15 
h^^ Mpc comoving for particle (rootgridsize) numbers 50'^, 
fOO^ 200^ The results are shown in Figure |H] We ran the 
Enzo simulations for 10 time steps starting ai z = 3.5. For 
Gadget-2 we ran the simulation from z = 3.5 to the same 
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redshift that Enzo had reached. The Enzo(AMR) simu- 
lations are about a factor 1.5-2 faster than the Gadget- 
2 simulations without star formation. Turning the AMR off 
leads to a speed-up of a factor five for the Enzo simulations. 
The Gadget-2 simulations without star formation spends 
most of its time calculating the hydrodynamics of the very 
high density gas, which for the Lyman-a forest is not neces- 
sary. Turning on the star formation in Gadget-2 leads to a 
spee d-up by a factor 30. As we have demonstrated here (see 
also lViel et alj (|2004) ) the effect of turning on star forma- 
tion in Gadget-2 on the statistics of the flux distribution 
is very small. 

Cosmos is an SGI Ahix 3700 with 152 Itanium2 (Madi- 
son) processors and 152GB of globally shared main memory. 
Each Madison processor has a clock speed of 1.3GHz, a 3MB 
L3 cache and a peak performance of 5.2 Gflops. The system 
is built from 76 dual-processor nodes, each with 2GB of local 
shared memory. These are linked by the SGI NUMAflexIII 
interconnects, which provides a high speed (3.2GB/sec bi- 
directional), low latency (sub- microsecond) network with a 
dual-plane, fat tree topology, connecting all processors with 
each other and with a single, globally shared and cache- 
coherent 152GB memory subsystem. We have only tested 
the two faster versions of the codes on the shared memory 
machine. The results are shown in Figure|H| Note the reversal 
in relative speed between Gadget-2 simulations with star 
formation and Enzo static grid simulation between the two 
architectures. Obviously Enzo ben efits more strongly f rom 
the shared memory architecture (see lO'Shea et all feOOSTl for 
a similar result). The Enzo(static) simulation run about 
a factor three faster on the shared memory machine than 
the Gadget-2 simulations with star formation. Reducing 
the refinement level in the AMR simulation may thus offer 
a good compromise between accuracy and speed for Lyman- 
a forest simulations with Enzo . 



4 DISCUSSION AND CONCLUSIONS 

We have performed a detailed comparison of Lyman-a for- 
est simulations with Gadget-2 , a TreePM-SPH code, and 
Enzo a Eulerian AMR code in order to asses the numeri- 
cal uncertainties due to a particular numerical implementa- 
tion of solving the hydrodynamical equations. The codes are 
similar with respect to the way in which they compute the 
gravitational forces at large scales but differ in the way they 
calculate gravitational forces on small scales; the codes use 
a Tree-PM and PM N-body algorithm, respectively. Their 
main differences lie, however, in the way in which they solve 
the gas hydrodynamics. Gadget-2 discretises mass using 
SPH methods while Enzo discretises space using adaptive 
meshes. The main results are as follows. 

• The differences in the dark matter power spectrum be- 
tween simulations with Enzo and Gadget-2 on scales rel- 
evant for measurements of the matter power spectrum from 
Lyman-a forest data are about 2% for an appropriate choice 
of box size and resolution. 

• The temperature density relation of simulations with 
Enzo and Gadget-2 differ very little. The PDF of the vol- 
ume weighted temperature differ by ~ 10 % probably mainly 
due to differences in the PDF of the gas density which are 
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Figure 8. A comparison of the CPU time required for the sim- 
ulation to run through a fixed redshift interval. Simulation pa- 
rameters are as described in the text and annotated on the plot. 
The thin lines are for a SUN distributed memory cluster and the 
thick linos are for a SGI shared memory (SM) machine. Only rel- 
ative values for the same architecture should bo considered. Note 
the reversal in relative speed between Gadget-2 simulations with 
star formation and Enzo static grid simulation between the two 
architectures. The thick dashed lino shows a linear scaling of CPU 
time with the number of particles for comparison. 

of the same order and at least partially caused by a slight 
mismatch in resolution. 

• The PDF of the flux distribution of simulations with 
Enzo and Gadget-2 agree very well. Typical differences are 
~ 5 — 10% probably again mainly due to a slight mismatch 
of the resolution. 

• The differences of the flux power spectrum of simula- 
tions with Enzo and Gadget-2 on scales relevant for mea- 
surements of the matter power spectrum from Lyman-a for- 
est data are about 5% for an appropriate choice of box size 
and resolution and simulations which fully resolve the Jeans 
mass. For simulations of lower resolution but larger boxsize 
the difference increase up to ~ 10%. Note that the differ- 
ences are scale and redshift dependent. 

Overall the Lyman-a forest simulations with Enzo and 
Gadget-2 agree astonishingly well. The choice of method 
for solving the hydrodynamical simulations appears to af- 
fect the gas distribution and its thermal state very little. 
It is also reassuring that two different implementations for 
solving the gravitational equations agree well. The corre- 
sponding uncertainties should contribute to the overall error 
budget of measurements of the matter power spectrum from 
Lyman-a forest data at the level of 3%. The total error in 
current measurement is signiflcantly larger and they should 
thus not be important. The main numerical uncertainties 
are instead due to a lack of sufficient dynamic range which 
typically makes correction of 5% for boxsize and resolution 
necessary. This will obviously improve as computational re- 
sources become more powerful. In practical terms memory 
requirements of simulations with Enzo without AMR and 
Gadget-2 are simUar. Enzo without AMR offers the high- 
est speed but requires somewhat larger corrections. Our re- 
sults suggest that if sufficient computational resources are 
available and sufficient care is employed the accuracy of nu- 
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merical simulations should not yet be a limiting factor in im- 
proving the accuracy of measurements of the matter power 
spectrum from Lyman-a forest data. 
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